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CALCULATION OF THE CENTERED ONE -DIMENSIONAL 
UNSTEADY EXPANSION OF A REACTING GAS MIXTURE SUBJECT TO 
VIBRATIONAL AND CHEMICAL NONEQUILIBRIUM 

By Laurence N. Connor, Jr. 

Langley Research Center 

SUMMARY 

A method of calculation based on a previously developed method-of-characteristics 
approach is presented for use in analyzing the nonequilibrium one-dimensional unsteady 
expansion of a reacting mixture of gases. The characteristic equations are written in a 
general form which permits the consideration of a number of rate processes. A multi- 
component gas model with a number of simultaneous rate processes is used, and both 
chemical and vibrational nonequilibrium are permitted. A procedure which utilizes the 
method of characteristics in a Lagrangian frame of reference is programed to yield solu- 
tions on the IBM 7094 electronic data processing system. Calculations are presented for 
a typical unsteady expansion to demonstrate the use of the program. 

INTRODUCTION 

The one -dimensional unsteady expansion has recently assumed new importance in 
the operation cycle of facilities designed to produce high-enthalpy flow for reentry simu- 
lation and chemical kinetic studies. The expansion tube discussed in references 1 and 2 
is designed to produce high-velocity flow for reentry simulation. In the operating cycle 
of the expansion tube an unsteady expansion is used to accelerate the test gas. This 
expansion is a critical phase in the expansion-tube cycle since any deviation from equi- 
libium at this point will alter the final state of the test gas. The need for a solution which 
is capable of completely evaluating the unsteady expansion is therefore apparent. 

For high-enthalpy gas flow, the role of chemical kinetics assumes new importance. 

In order to analyze such flow accurately, a knowledge of the rates of the reactions involved 
is necessary. The unsteady expansion can provide a means of processing a gas to estab- 
lish conditions under which chemical recombination rates can be directly measured. Such 
an application of the unsteady expansion is discussed in references 3, 4, and 5. The oper- 
ating conditions involve temperatures and pressures at which internal energy excitation 
and dissociation of the molecular species of the processed gas occur. The rate at which 
these processes occur appreciably affects the flow properties of the gas processed by the 



expansion. An analysis of the unsteady expansion in which rate -controlled phenomena can 
be included and for which the effects of such processes can be ascertained is therefore 
required. 

The first approach to predicting the structure of a nonequilibrium unsteady expan- 
sion is found in reference 6. A centered rarefaction wave in which only vibrational 
nonequilibrium is present is treated. That work was extended in references 7 and 8 by 
introducing the "ideal dissociating gas" of references 9 and 10 into the treatment of 
dissociative nonequilibrium. The calculation procedure and computed results for the 
unsteady expansion of a vibrationally relaxing nitrogen-oxygen mixture are presented 
in reference 11. The numerical approach, originated in reference 6 and used in refer- 
ences 7 and 11, is utilized in the present work. 

The basic objective of the present investigation is to provide an analysis of the one- 
dimensional unsteady expansion of a reacting mixture of gases subject to vibrational and 
dissociative nonequilibrium. The physical model of the expansion is the constant -velocity 
withdrawal of a piston from a cylinder of gas which is initially in chemical equilibrium. 
The chemistry and thermodynamics of the problem are formulated generally so that vari- 
ous gas mixtures and reactions may be considered. The equations of motion governing an 
unsteady expansion are used to establish characteristic lines and compatibility relations 
along the lines. 

A procedure is established for calculation of the flow field by the method of charac- 
teristics and is them programed in FORTRAN language to yield solutions on the IBM 7094 
electronic data processing system. The solution of a typical unsteady expansion in an 
oxygen-nitrogen mixture is presented to illustrate the general structure of such an expan- 
sion. The variations of thermodynamic properties, velocities, and compositions through 
the expansion fan are depicted for gas particles having different residence times in the 
expansion, and thus the effects of nonequilibrium on the expansion are illustrated. 

SYMBOLS 


Al,A2,A 3 

Ai,Bi,Ci 

a 

b 

D j’ E j’ F j 


dimensionless parameters defined by equation (18) 
constants used in vibrational-relaxation-time expressions 
speed of sound 

Lagrangian coordinate denoting a fixed mass element 
chemical rate constants 
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internal energy per unit mass 


factor used to provide exclusion and inclusion of rotational energy of 
monatomic and diatomic molecules, respectively (fi = 0 for monatomic 
molecule; fi = 1 for diatomic molecule) 

enthalpy per unit mass 

species identification 

reaction identification 

total number of rate -controlled variables 

rate-controlled-variable identification 

specific reaction rate coefficient for jth reaction 

catalytic species in chemical reaction 

total number of reactions 

total number of species 

Avogadro number 

pressure 

rate-controlled variable considered in appendix A 
general rate-controlled variable (k = 1, . . .,K) 
universal gas constant 
translational temperature 
time 


reference time 



U piston velocity 

u velocity 

x position coordinate of fluid particle 

Xj mass fraction of ith species for i = 1, . . .,n (xq represents the mass 

fraction of atomic oxygen) 

y modified Lagrangian coordinate, b/t 

y frozen specific heat ratio of the undisturbed gas 

heat of formation of ith species, energy per molecule 

6i characteristic vibrational temperature of ith species 

H molecular weight of mixture 


Mi 


molecular weight of ith species 


'i] 


stoichiometric coefficient of ith species of reactant in jth reaction 


i] 


stoichiometric coefficient of ith species of product in jth reaction 


density 


vibrational energy per mole of ith species for i = 1, . . .,n anc * CT C>2 

represent vibrational energies of N 2 and O 2 , respectively) 


vibrational relaxation time of ith species 


u> 


rate of change of rate-controlled variable q 


u> k 


rate of change of kth rate-controlled variable 


Subscripts: 

A,B,C,D,P points in calculation procedure defined in figure 9 
AP,BP,DP denote average of values at two grid points 
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c required for complete frozen expansion to p = 0 

eq equilibrium conditions 

f results for frozen flow 

n evaluated at intermediate condition (fig. 1) 

o initial condition in undisturbed gas 

p evaluated with pressure held constant 

q evaluated with rate -controlled variables in appendix A held constant 

q^ evaluated with rate -controlled variable held constant 

p evaluated with density held constant 

+,-,0 characteristic directions 

A bar over a symbol denotes a nondimensional quantity. 

THEORETICAL DEVELOPMENT 
Physical Model of Expansion 

A physical picture of a centered one-dimensional unsteady expansion can be obtained 
by considering figure 1. For illustrative purposes the structure of the expansion is shown 
for a perfect gas. A semi-infinite reservior of constant cross-sectional area contains a 
gas mixture assumed to be in a state of chemical and thermodynamic equilibrium. At 
time zero, the piston is withdrawn with constant velocity, and a rarefaction wave which 
propagates into the gas in the reservoir is produced. 

The expansion can be divided into two regions. The first is referred to as the 
expansion fan. In this region, temperatures and pressures decrease rapidly during the 
passage of the rarefaction wave. If the reservoir gas exists initially in a dissociated 
state, the atoms must recombine as they are processed by the expansion if chemical 
equilibrium is to be maintained. In like manner, the vibrational energy must decrease. 
Both of these phenomena are controlled by rate processes. 

The effects of rate processes can be understood by considering the different particle 
paths in figure 1. Path I represents a particle which is very near the piston face initially. 
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Its residence time in the expansion fan is very short and the rate -controlled quantities 
have insufficient time for adjustment to equilibrium. Since very little adjustment occurs, 
the particle essentially experiences a frozen expansion. Particle HI has a relatively long 
residence time in the expansion and the rate processes are permitted adequate time for 
the flow to relax to equilibrium. Such a particle undergoes an equilibrium expansion. 

The energy contained in dissociation is returned to the translational and rotational energy 
modes of the gas. Particle path II represents nonequilibrium flow, an intermediate region 
between the two extremes. The residence time is of a magnitude that allows some relax- 
ation of the rate -controlled quantities toward equilibrium. The present investigation is 
directed toward the analysis of this nonequilibrium region. 

Governing Equations 

The Lagrangian coordinate system can be used effectively in the present problem. 
Such a system has the advantage of readily establishing a particle path along which the 
rate and energy equations can be applied. In the Lagrangian frame of reference, attention 
is given to what happens to each individual fluid particle in the course of time. Each fluid 
particle must therefore be labeled. This labeling is accomplished by designating each 
particle according to its position coordinate x in the x-t plane at some reference time 
t = t 0 . This value is called b, the Lagrangian coordinate specifying the particle being 
considered. The two independent variables in the Lagrangian system are therefore b 
and t. 

Flow equations .- The equations expressing the conservation of mass, momentum, 
and energy for one-dimensional unsteady expansion in the Lagrangian coordinate system 
are presented in references 12 and 6 and are as follows: 


9b p 2 9t 

(1) 

+ J_M = 0 

9t P 0 9b 

(2) 

9h 1 9p _ 0 
8t ” P 9t 

( 3 ) 


Equation of state .- The state of the reacting gas mixture can be expressed in the 
general form 

h = h(p,p,qj,. . .,q K ) (4) 


6 


The variables denoted by for k = 1,. . .,K are used to specify the vibrational 
energy and/or the chemical state of the gas. In actuality, these variables express the 
vibrational energy and/or the mass fraction of each component species in the gas mix- 
ture, that is, 


qi = q£ = o 2 , 


•> q n = CT n 


‘In+l _ x l> 1n+2 _ x 2» • • •> *l2n ~ x n 


( 5 ) 


where is the vibrational energy per mole of the ith species for i = 1,. . ,,n, xj is 
the mass fraction of the ith species for i = 1,. . .,n, n is the total number of species, 
and K is the number of nonequilibrium variables (K = 2n). For convenience in handling 
the equations in the present paper q^ is used to represent the vibrational and dissocia- 
tive state of the gas. When more detailed operations are required, equation (5) is intro- 
duced so that each variable can be treated in the proper manner. 

Rate equations . - For the nonequilibrium expansion, the value of q^ at any point in 
the flow depends on the rate equations and the process which the gas has undergone. The 
rate equations can be expressed in the general form 


where 



w k = “>k(p»p»qi»- • -^k) 


( 6 ) 


and are subsequently presented in detail for each rate process considered. 


Characteristic Equations 

The characteristic equations in the Lagrangian coordinate system are derived in 
appendix A. Along the + and - characteristic lines 





the respective compatibility relations are 
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I 



d£ 

dt 


d£ 

dt 


K 


9h 


X 8q. Wk 


+ paf 


du o k— 1 k 

3r“ 


dt 


a f 


ah 

ap 


K 

z 

du o k=l 
' P f dT " af 


8h 

0< lk 


w k 


ah 

ap 




The particle path in the Lagrangian system is given by 


( 8 ) 



and, along this path, both the energy equation 


and the rate equation 


dh 1 dp _ 0 

dt P dt 



0) 


( 10 ) 


(ID 


must be satisfied. 

Writing the characteristic equations in dimensionless form facilitates the computa- 
tion necessary for a numerical solution. Thus, the following dimensionless quantities are 
defined: 


af 

af = 

a f,o 


u- u 


U r 




p = -E_ 

P Po 


a f,ot' 


h = f 

h 0 




( 12 ) 


The subscript o refers to the value of the quantity in the undisturbed gas before the 
piston is withdrawn, and U c denotes the piston velocity required for a complete frozen 
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expansion to p = 0. The reference time t' can be chosen as necessitated by the par- 
ticular expansion being calculated, with the provision that it not be zero. 


These expressions are used to write the dimensionless characteristic equations of 
an unsteady expansion in the Lagrangian coordinate system. Along the + and - char- 
acteristic curves 



the ordinary differential equations 


<fc + Al= L-d£>- A 2 f 
dt p af dt 


Al J r * + A 2 f 

dt paf dt 


- K _ ^ 

4y^ k =o 

p<* 4 9q * 

9p 

K 

a f Y 9h n 

9 p > 




(13) 


(14) 


must be satisfied, and along the particle path 



(15) 


the energy equation and the rate equations may be written as 


dh A 3 dp _ Q 
dt P dt 


and 


where 



Al = 


Po 


P o%o V c 


A 2 = 



U c 


a 3 = 


Po 

^qPq 


(17) 


(18) 


9 



The application of the boundary conditions to the problem can be implemented by a 
modification of the coordinate system. This modification, first used in reference 6, is 
accomplished by dividing the Lagrangian particle coordinate b by the time variable t 
to define a new variable y — that is, y is defined by the relation 

y = * (19) 

t 

The independent variables in the modified Lagrangian system are therefore y and t. 

Use of such a coordinate system eliminates the singularity which exists as t — 0 and 
allows the boundary conditions to be applied. In addition, the network of grid points 
becomes more nearly rectangular in the y,t coordinates. 

When b is replaced by yt in equations (13) and (15), the compatibility equations 
remain the same and the equations for the characteristic curves become 


= P»f 

and 

= 0 

0 




( 20 ) 


( 21 ) 


Equations (20), (14), (21), (16), and (17) make up the set of characteristic equations 
which, after conversion to finite-difference form, are used to obtain a numerical solution. 


Boundary Conditions 

Three conditions must be prescribed to determine the solution. Of these conditions, 
two are used in the calculation of the expansion fan and the other is used in the analysis of 
the near -steady region between the tail of the expansion fan and the piston face. The 
physical reasoning supporting the selection of the prescribed conditions is presented. 

Two initial conditions are prescribed for the expansion region. First, the undis- 
turbed gas in the reservoir is assumed to be in a state of thermodynamic equilibrium. 

The application of this condition in each coordinate system can be understood by consid- 
ering figure 2. In y,t coordinates, this first boundary condition can be written in the 
form 
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( 22 ) 



The second boundary condition is obtained by considering the expansion as t — 0. 

A particle processed by the expansion for this second condition has a residence time 
which also approaches zero. The particles adjacent to the piston face experience such an 
expansion but, because of the short residence time, no adjustment of rate-controlled 
quantities occurs. It is shown mathematically in references 3 and 6 that, as t 0, the 
properties throughout the process are those of a frozen expansion and can be expressed 
as a function of y. In the dimensionless modified Lagrangian coordinate system, this 
boundary condition may be written as 

2y 0 

p(y,o) = p f (y) = y r ° +1 
2 

p(y,o) = P f ( y) = y y ° +1 > (23) 

2Vi 

u(y,o) = u f (y) = u Q + 1 - y r ° +1 

<*k(y’°) = %,f 

The values of £, of course, remain constant at the initial undisturbed values of q^ 
for the frozen boundary conditions. 

A final boundary condition is needed to permit the calculation of the region of near- 
steady flow between the expansion fan and the piston face. The condition required is that 
the gas velocity at the piston face be the same as that of the piston. Such a condition 
must exist if no cavitation occurs. Cavitation, described as a breakdown in continuum 
flow, occurs only if the piston velocity is greater than the limiting velocity that the gas 
could attain in an expansion to zero pressure. Thus, by choosing the piston velocity to be 
less than U c , the possibility of cavitation is eliminated. 

Thermodynamic Model of Reacting Gas Mixture 

The equation of state and the rate equations have been taken to be of the general 
forms given by equations (4) and (6), respectively. It is necessary next to develop spe- 
cific relationships for the thermodynamic properties of the reacting gas and to establish 
detailed rate expressions for the various nonequilibrium processes. 
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Thermal equation of state.- The reacting gas is assumed to be made up of a mix- 
ture of perfect gases. The pressure of the mixture is equal to the sum of the partial 
pressures of the component species. The thermal equation of state can therefore be 
written as 


P = 


pRT 


where 



( 24 ) 


( 25 ) 


and where R is the universal gas constant and x^ and are the mass fraction and 
molecular weight, respectively, of the ith species. 

Internal energy.- The internal energy of the gas contains the contributions from 
translational, rotational, vibrational, and electronic energy modes plus the energies of 
dissociation and ionization. Since a mixture of perfect gases is considered, the internal 
energy of the mixture is equal to the sum of the internal energies of the component 
species. The internal energy per unit mass e can be expressed as the sum from 
i = 1,. . .,n of the translational, rotational, vibrational, and electronic energies and the 
energies of dissociation and ionization. 

A detailed evaluation of the contributions of the various energy modes can be made 
by using the appropriate partition functions. Several assumptions are made in the formu- 
lation of the expression for internal energy. The gas mixture is restricted to monatomic 
and diatomic species. Electronic excitation and ionization are neglected. Compositions 
and vibrational energies are assumed to be rate controlled. By substituting the various 
contributions to the internal energy into the summation equation for e, the internal 
energy of the reacting mixture of gases is seen to be 


n 

6 = I S| RT + fi ( RT + CT i) + N ° A i 
1=1 lL 


( 26 ) 


where T is the translational temperature of the gas, N 0 is the Avogadro number, and 
Aj is the heat of formation in energy per molecule for the ith specie. The factor fi is 
used to permit the existence of rotational and vibrational energies for molecular species 
and to eliminate these energy modes for atomic species. For diatomic molecules, 
fi = 1; for monatomic species, fj = 0. 
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Enthalpy. - By using the definition of enthalpy 


h = e + E 
P 

and equations (24), (25), and (26), the expression for the enthalpy in units of energy per 
unit mass of the mixture becomes 


h = 



( 27 ) 


The assumption was made previously in this analysis that the caloric equation of state 
could be written as equation (4) - that is, in the form 


h = h(p,p,q v . . .,q K ^ 


The variable q k was used to represent all the rate-controlled variables - in this 
instance, the mass fraction Xf and the vibrational energy of of each species of gas. 
Equation (27) expresses the enthalpy in the prescribed form and permits the evaluation 
of the enthalpy derivatives which appear in the characteristic equations. 


Enthalpy derivatives .- Certain enthalpy derivatives which appear in the charac- 
teristic equations must be evaluated in terms of the gas model. The frozen speed of 
sound SLf, which appears in both the characteristic and the compatibility equations, is 
expressed in reference 13 as 


a f 



(28) 


In order to 

p, p, and 
are seen to 


8h 


evaluate af it is necessary to express — and 

^ apy P,Qk 

q^. By using the enthalpy expression given in equation 
be 



| in terms of 
Ik 

these derivatives 





(29) 
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The frozen sound speed can then be written as 



( 30 ) 


It can readily be shown that the quantity in the brackets is the ratio of the frozen specific 
heat defined in the conventional way. 

K 

In the compatibility equations, the term / must be evaluated. Since 

k=i Sqk 


K n n 

v = yjto<&i + v_9hdci 

L 9q, k L 0x, dt L am dt 

k=l k i=l 1 i=l 1 


( 31 ) 


expressions for the derivatives — and — are needed. For the present gas model, 


00* 


these expressions are found to be 


_0h 

0x 


L = (| +f w + f a./| + fi ) 

i V 2 JPH H H H H\ 2 7 


( 32 ) 


and 


x £i 

'i ,X i 


8h _ 
BOk 


( 33 ) 


K 


When dissociative and vibrational nonequilibrium exists, the term ) — — may then 

i—i dq. 

k=l K 

be expressed in the following form: 

V 8h Y f?5 + f \PM + f i CT i + N oAi Y X i/ S + f.^l dX i + f X i fidai r 34 ) 

Z % Wk= Z (2 + f te + + — -p-p Z MTl2 + f Vdr + Z df (34) 

k=l i=l [j i=l J i=l 

dxs do* 

The rate equations, discussed subsequently, provide expressions for and 

in terms of thermodynamic properties, compositions, and vibrational energies. 
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Rate eq uations . - The rate equations were postulated to take the form of equa- 
tion (6) — that is, 



This expression must now be expanded to more specific terms. 

A rate equation is first established to describe the dissociation and recombination 
rates of the reacting gas species. By using the rate expression presented in reference 14 
in terms of the mass fraction, the rate of change of the mass fraction of the ith species 
due to all reactions may be expressed as 


m 



j=l 


cbq 

dt 


m 


i 


j=i 




(35) 


where x^ is the mass fraction of the ith species, j denotes the reaction, is the 

stoichiometric coefficient of ith species of reactant in the jth reaction, is the 

stoichiometric coefficient of ith species of product in the jth reaction, and kj is the 
specific reaction rate coefficient for the jth reaction. The specific reaction rate coeffi- 
cient is assumed to be a function of temperature only and is an empirically determined 
coefficient which can be written in the following form: 


k j 



(36) 


The form of the chemical rate equation given by equation (35) is adequate for both forward 
and reverse reactions, provided each is treated as a separate reaction. 

The adjustment of the vibrational energies toward equilibrium is also controlled by 
a rate process. In the present investigation, coupling between dissociation and vibra- 
tional energies is not considered. In addition, any molecule formed by the recombination 
process is assumed to possess the same average vibrational energy as other molecules 
in the flow. The vibrational rate equation may then be written as 

dOj CTi prt — Oi 

i = i (37) 

dt 

where the vibrational relaxation time 


takes the form 
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(38) 


AjT Bi exp I 


T i = — 


vT 1 / 3 / 


1 - expl-^ 


which readily fits available empirical data and agrees with the form theoretically pre- 
dicted in reference 15. The equilibrium vibrational energy is assumed to be given by the 
expression 


R0i 



(39) 


Vibrational-equilibrium option .- In some situations the vibrational relaxation is 
much more rapid than the recombination processes. It is therefore advantageous to 
assume vibrational equilibrium and treat the dissociation and recombination reactions as 
the only rate processes. For conditions considered in this investigation the variation of 
equilibrium vibrational energy with temperature can be represented by a series of 
straight lines with no appreciable loss in accuracy. Such an approximation simplifies 
the evaluation of enthalpy derivatives. Assuming the vibrational energy to be essentially 
in equilibrium permits more rapid calculation than could be attained by considering vibra- 
tional energies as rate controlled. 


SOLUTION PROCEDURE 


In appendix B the detailed operation of converting the nondimensionalized character- 
istic equations to finite -difference form and of establishing a procedure for calculating the 
desired flow region is presented. The calculation procedure described in appendix B was 
programed in FORTRAN IV language for computation on the IBM 7094 electronic data 
processing system. Capabilities and limitations of the program are discussed herein. 

The constituents of the gas mixture must be monotomic or diatomic species if the 
present thermodynamic model is to be used. Conditions throughout the expansion must be 
such that ionization and electronic excitation can be neglected. 

The program is designed to handle a number of chemical species and the associated 
chemical and vibrational rate equations. Computer time and storage limitations make 
difficult the treatment of simultaneous rate processes in which the overall rates of the 
important relaxation processes differ by several orders of magnitude. Since the grid 
size must be small enough to accommodate the fastest rate, the number of grid points 
required for accurate characteristic calculation places a practical limit on such work. At 
the present time no suitable method, comparable to the one introduced in reference 16 for 
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use in one -dimensional flow, has been developed for characteristic calculations. This 
grid-size limitation, coupled with the assumptions incorporated in the thermodynamic 
model, restricts the use of the computational program to fairly simple gas mixtures. 

These restrictions do not preclude the utilization of the program in many practical 
applications. 

The program can be used to calculate the flow field over a large portion of the non- 
equilibrium region between the limiting cases of frozen and equilibrium flow. For flows 
expanded to very low densities, the number of grid points necessary to calculate the entire 
region presents a problem because of the excessive computing time required. The flow 
in the initial portion of the expansion reaches equilibrium rapidly; however, as the expan- 
sion proceeds, the density drops with a resulting decrease in the recombination rate. For 
expansions to low densities, a large region must therefore be calculated in the t direc- 
tion of the y-t plane to follow the flow to equilibrium. Some computational problems 
have been encountered in the near -equilibrium region. The procedure used in the charac- 
teristic calculation provides a natural perturbation of the nonequilibrium variables from 
equilibrium at the beginning of the calculation but, as calculation proceeds in the t direc- 
tion and equilibrium is approached, numerical inaccuracies may occur because of the 
nature of the rate equations. The approach in the present investigation has been to com- 
pare the compositions and vibrational energies with the equilibrium values at the local 
conditions and, if they are sufficiently close, to equate all further values in the t direc- 
tion to the equilibrium values. This approach has proved adequate for the simplified gas 
model discussed previously. 

The program might be further improved by using a method similar to that pre- 
sented in reference 17 to aid in the calculation of the near -equilibrium region. It might 
also be noted that, for many useful applications of the program, the calculations do not 
have to be carried to the equilibrium limit. 

The program may be used for expansions in which the pressure drop does not 
exceed two orders of magnitude. For most applications of current interest, the rate- 
controlled quantities will have essentially frozen or become insignificant for this degree 
of expansion and any further expansion can be approximated by frozen flow. 

Results calculated for the region of expanding flow may be printed out in two differ- 
ent ways. The thermodynamic properties, velocities, vibrational energies, and mass 
fractions can be printed either at each grid point in the y-t plane or along specified 
particle paths. A numerical integration routine is also included to calculate the distance 
traveled by the specified particles. 



RESULTS AND DISCUSSION 


Results of Sample Calculation 

Calculation of the unsteady expansion of a gas experiencing chemical and vibrational 
nonequilibrium was made for a specific set of initial conditions. Results are presented to 
demonstrate the utility of the computational program and to point out the salient features 
of such an expansion. 

A simplified air model, made up of N 2 , O 2 , and O, is used in the analysis. The 
oxygen dissociation-recombination reaction 

O 2 + M r 20 + M 

where M is O, 0 2 , or N 2 , is assumed to control the change in chemical Composition- 
Different rate constants may be used for each catalytic species M. The chemical rate 
constants required in equation (36) were obtained from reference 18 and are listed in 
table I. 

Rate equations of the form given by equation (37) are used to describe the relaxation 
of the vibrational energies of molecular oxygen and nitrogen. The vibrational-relaxation- 
time constants were obtained by curve fitting the values of reference 19 and are given in 
table n. Recent experiments described in references 20 and 21 indicate that vibrational 
relaxation rates in expanding N 2 and air can be expected to be 15 times as fast as those 
observed in compressive flows. Oxygen vibrational deexcitation, however, was measured 
in an unsteady expansion (ref. 5) and the results were found to agree with calculations 
made by using rates measured in a compressive flow. Some question exists, therefore, 
as to the correct rate data to use in the present investigation. Although the data of refer- 
ence 19 predict rates which are possibly too slow in this expanding-flow situation, they 
are used herein to illustrate the handling of simultaneous vibrational and chemical rate 
processes. Use of the faster rates would result in vibrational equilibrium throughout the 
expansion. No exchange of vibrational energies between N 2 and O 2 molecules, as dis- 
cussed for compressive flow in reference 22, is considered in the present study. The 
further assumption is made that no coupling exists between the chemical and vibrational 
rate processes. 

The undisturbed gas, before being processed by the expansion, was assumed to 
exist in a state of chemical equilibrium at a temperature of 3500° K and a pressure of 
1.0968 X 10® dynes/ cm2. jp or suc h conditions, the mass fraction of NO is small, the 
amount of N is negligible, and the other assumptions used in the formulation of the gas 
model are met. The simplified air model therefore appears adequate for the present 
application. 


18 



Variations in compositions, vibrational energies, and flow properties along speci- 
fied particle paths through the unsteady expansion are presented in figures 3 to 8. The 
reduced time coordinate t/t Q , used as the abscissa in these figures, is obtained by 
dividing the actual time the particle spent in the expansion by the time at which that par- 
ticle would have entered a frozen expansion. Property values along different particle 
paths for a given value of t/to describe particles which have different residence times 
in the expansion but which have been processed by approximately the same fraction of the 
total expansion. Specific particle paths are designated by t Q , the time the particle enters 
the fan. This value is also a direct measurement of the initial distance of the particle 
from the center of the expansion at time zero. Equilibrium properties are calculated by 
the procedure described in reference 23. 

Composition changes .- The nature of the variation in composition, as depicted in 
figure 3, is considered. As the value of t Q of the particles increases (an indication of 
increased residence time in the expansion), the amount of recombination increases and 
thus the mass fraction of atomic oxygen in the flow decreases. The atomic oxygen 
decreases very rapidly during the initial part of the expansion and then levels off as if it 
were asymptotically approaching a constant value. This phenomenon is similar to the 
composition freezing which occurs in hypersonic nozzles and can be attributed basically 
to two factors. First, a decrease in oxygen atoms decreases the recombination rates. 

As recombination occurs, the mass fraction of atomic oxygen in the flow decreases. The 
recombination rates, which are proportional to the square or cube of the mass fraction of 
atomic oxygen, therefore decrease rapidly. The second factor influencing the recombina- 
tion rates is the density. As expansion occurs, the density drops rapidly. From the 
chemical rate equation, it is apparent that the recombination rate is proportional to the 
square of the density. Thus, a density decrease produces a rapid decrease in the rate of 
recombination. This rate decrease due to density change more than offsets the effect 
produced by the inverse temperature dependence of the reaction rate coefficients. The 
net result of the property changes produced by the expanding flow on the recombination 
rates is a sizable reduction in these rates. 

Vibrational deexcit ation.- The variations in the vibrational energies of molecular 
nitrogen and of molecular oxygen are displayed in figures 4 and 5, respectively. Again a 
tendency toward freezing is noted as the gases expand to low densities. Of particular 
interest is the behavior of the oxygen vibrational energy. The initial rate of deexcitation 
is faster for the particle with the smaller residence time. The flow conditions experi- 
enced by this particle are far from equilibrium with respect to dissociation and nitrogen 
vibrational energy. Very little chemical or vibrational energy has been returned to the 
flow. The flow, therefore, possesses a low translational temperature causing the value 
of the term cq e q - cq in the oxygen vibrational rate equation to be large and thereby 
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providing the necessary potential for producing rapid oxygen vibrational deexcitation. As 
the expansion proceeds, the pressure-dependent relaxation time begins to dominate the 
oxygen vibrational energy relaxation, and the tendency toward freezing begins. The par- 
ticle more distant from the corner has more residence time for relaxation and thus its 
vibrational energies approach the frozen state more slowly. 

Thermodynamic pr operties.- The translational temperature is the thermodynamic 
property influenced to the greatest extent by the existence of nonequilibrium. Figure 6 
shows the temperature profiles along various particle paths through the unsteady expan- 
sion. As recombination and deexcitation of vibrational energy modes occur, energy is 
returned to the flow and manifests itself partially in the form of a rise in translational 
temperature. Temperature is therefore seen to be highly dependent upon the chemical 
and vibrational state of the gas. The pressure variation along specified particle paths is 
displayed in figure 7. For a high degree of expansion, the nonequilibrium static pressure 
can deviate appreciably from the equilibrium value and can serve as an indicator of the 
degree of nonequilibrium of the flow. Density variation is not shown graphically because 
of its slight sensitivity to the chemical and vibrational state of the gas. 

Velocity variation . - The flow velocity is of considerable interest in many applica- 
tions of the unsteady expansion, particularly in those dealing with the design of facilities 
for producing high-velocity, high-enthalpy flow. The effect of the chemical state of the 
gas on the velocity is therefore very important. It can be seen from figure 8 that, as the 
recombination increases, the velocity of the gas increases. This increase in velocity is 
to be expected inasmuch as some of the energy released by the recombination is converted 
into flow energy. In order for the maximum velocity to be obtained in an unsteady expan- 
sion, the process must take place in a state of chemical equilibrium. 

Proposed Applications of Program 

The practical flow situations which can be calculated by using the present computa- 
tional program can be divided into two basic categories. The sample calculation just pre- 
sented provides an example of the first category. This category of application deals with 
the determination of the effects of nonequilibrium on the unsteady expansion of a gas being 
processed for use as a test medium. The objectives of the analysis are to predict flow- 
property variations through the expansion and to ascertain the state of the test gas. Tem- 
perature and pressure ranges must be restricted to conform with the assumptions of the 
thermodynamic model. The number of species and reactions considered must be held to a 
minimum to make computation practical. For many problems of this nature, such as the 
analysis of expansion-tube flow (ref. 1), these restrictions do not prohibit realistic 
calculation. 
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The second category in which the present program can be utilized is the analysis of 
unsteady expansions in kinetic rate studies. The use of the unsteady expansion in pro- 
ducing conditions for which recombination and vibrational relaxation can be measured is 
described in references 3, 4, and 5. An experiment of this kind would be performed with 
a mixture containing a large percent of an inert gas such as argon and a smaller amount 
of a gas such as nitrogen or oxygen for which the rate is to be measured. The present 
program is well suited for calculating such a flow. The calculations would be used to 
predict the ideal test conditions for the experiment and would provide the basis for com- 
parison between measured and calculated values needed to determine actual values of rate 
constants. 


CONCLUDING REMARKS 

A procedure, based on the method of characteristics applied in a Lagrangian frame 
of reference, has been presented for calculating the centered one -dimensional unsteady 
expansion of a multicomponent gas mixture subject to chemical and vibrational nonequi- 
librium. A computer program was developed to perform the calculation on the IBM 7094 
electronic data processing system. The results from a sample calculation are presented 
to show the utility of the program and to illustrate the pertinent features of a nonequi- 
librium unsteady expansion. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., September 2, 1966, 

129-01-08-16-23. 
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APPENDIX A 


DERIVATION OF CHARACTERISTIC EQUATIONS 


Nonequilibrium Characteristics 

Equations (1), (2), (3), and (6), written in Lagrangian coordinates, describe the one- 
dimensional unsteady flow of a reacting gas mixture. This system of equations is hyper- 
bolic in nature and therefore has real valued characteristics. The equations defining the 
characteristic curves and compatibility relations are developed in this appendix by using 
these basic governing equations. In order to simplify the ensuing manipulations, only one 
of the rate -controlled variables q^ is considered and is denoted simply as q. It is 
apparent subsequently that this simplification does not limit the generality of the develop- 
ment of the characteristic equations. 

First, the caloric equation of state, equation (4), is used to eliminate derivatives of 
p. This equation may be written in the following form: 

ah = /9h\ 0 £ + /8h\ 8£ + /0h\ ^ ( A l) 

^ W) PA Bt W p , q « Wp, p W 


By introducing the energy equation, equation (3), and the frozen-sound-speed expression, 
equation (28), equation (Al) can be rearranged to give 
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Equation (1) can then be written as 
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Since p, u, h, and q are functions of the independent variables b and t, the 
following equations may be written: 
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Characteristic curves for the set of governing equations — that is, equations (2), (3), (6), 
(A3), (A4), (A5), (A6), and (A7) — can exist only if the determinant of the matrix of coeffi- 
cients of this set is equal to zero. The characteristic directions are obtained by equating 
the following determinant of coefficients to zero and evaluating it for that condition: 
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This equation reduces to 



It should be noted that, for each additional value of q considered, the factor ^ 

dt 

is raised to one higher power. An increase in the number of rate processes considered 
does not, therefore, introduce any additional characteristic directions. 

By setting each factor of equation (A9) equal to zero, the characteristic directions 
are seen to be 
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Equations (A10) and (All) are given in the text as equations (7) and (9). 

The compatibility relations which apply along these characteristic curves may now 
be easily obtained. Along the directions in equation (A10), equations (2) and (A3) can be 
combined to produce the following ordinary differential equations: 
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Equation (A12) is expressed in the text 

Along the particle path given by 
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as equation (8). 
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and 


dh 1 d£ = o 
dt " P dt 


(A13) 


^ = w(p,p,q) 


(A14) 


respectively. In the text, equation (A13) is given as equation (10) and equation (A14) is 
given in more general form as equation (11). 

It is interesting to compare the characteristic equations presented with similar 
equations for a perfect gas. The perfect-gas characteristics show that, along the lines 
given in equation (A10), the compatibility relations are 


(A15) 

du 1 dp _ Q 
dt pa dt 

In the nonequilibrium expansion the frozen speed of sound aj assumes the role of the 
perfect-gas unambiguous sound speed a, so the directions of the nonequilibrium and 
perfect-gas characteristic lines are similar. The compatibility equations, however, are 
not of identical form because of the addition in the nonequilibrium expansion of the term 
which is due to the inclusion of rate processes in the formulation of the problem. 

Frozen Characteristics 

In the limiting case of frozen flow, the rate -controlled quantities have no time in 
which to relax, q^ is constant, and therefore = 0. The frozen flow then reduces to 
the perfect-gas expansion in which the frozen compositions and vibrational energies are 
used in the calculation of the sound speed. 

Equilibrium Characteristics 

If the flow is to remain in thermodynamic equilibrium throughout the expansion, 
all rate processes must be assumed to occur infinitely fast. The characteristic equa- 
tions for the equilibrium expansion cannot be obtained as a limiting case of the nonequi- 
librium process since the rate equations become indeterminate. The characteristic 
equations can be developed separately for the mathematical model of equilibrium. The 
equilibrium characteristic equations are equivalent to the characteristic equations for a 
perfect gas with a perfect-gas sound speed replaced by the local equilibrium speed of 
sound a e q. 


du 

dt 
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NUMERICAL CALCULATION PROCEDURE 


The calculation procedure presented in this appendix closely follows the method 
outlined in reference 7. It consists of the conversion of the characteristic equations to 
finite -difference form and the establishment of a step-by-step computational scheme. 

Finite-Difference Form of Characteristic Equations 

When the nomenclature of figure 9 and second-order finite -difference approxima- 
tions are used, the equations for the characteristic curves given by equation (20) become 

yptp - Wa = Nap(*P - T a) c b1 ) 

and 

y P Fp - y B F B = -[p a f] Bp (Fp - F b ) (B2) 

and the respective compatibility relations given by equation (14) become 


up - u A + Ai 





and 


(B3) 


Up - u B 




(B4) 


A bracketed expression denotes the average of the value of the term contained in the 
bracket at the two grid points denoted by the bracket subscript. In finite-difference form 
the particle path, equation (21), becomes simply 


Yptp = 

Along the particle path, the energy and rate equations may be written as 


(B5) 
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and 
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Equations (Bl) to (B8) form the basis of the characteristic program for obtaining a 
numerical solution. 

Calculation Procedure for Expansion Fan 

For the expansion region, the calculations are carried out in the dimensionless 
modified Lagrangian coordinate system. The nature of the grid used in the calculations 
is shown in figure 10. 

From the initial conditions it is known that, along the line t = 0, all conditions are 
those of a frozen expansion. Thus, the values of all quantities at points (M,l) are known. 
Also, the properties along the line y = 1 are those of the undisturbed gas. Thus, all 
quantities are known for points (1,N). The starting point for the calculation is obviously 
point (2,2) inasmuch as points (2,1), (1,1), and (1,2) are known from the boundary condi- 
tions. The finite-difference equations are used to obtain the coordinates, properties, 
compositions, and velocities at point (2,2). Once point (2,2) has been calculated, 
points (2,2), (1,2), and (1,3) are used to calculate point (2,3). Thus, the calculation 
marches up the (2,N) line until the limiting N is reached. Then by using points (3,1), 
(2,1), and (2,2), point (3,2) is calculated and the same marching process is used along 
line (3,N). A marching process is therefore established which moves to the left after 
each successive column is calculated. This process continues until all desired points 
have been calculated. 

The preceding discussion outlines the general calculation of the entire flow field. 
The key operation is the calculation of values at any unknown point by using the known 
values at three adjacent points and the finite -difference characteristic equations. This 
operation is an iterative process and requires a fairly elaborate calculation procedure. 
The following steps are used in this calculation: 

1. All quantities are known at points A, B, and C (see fig. 9) from the boundary con- 
ditions or previous calculations. As a first assumption, the values of the thermodynamic 
properties, compositions, and velocity at point P are equated to the respective values at 
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point A. This assumption is accurate because, for flow with no nonequilibrium present, 
a simple flow region would exist and properties at points A and P would be identical. Dif- 
ferences in properties at points A and P are due then to nonequilibrium effects and not to 
the expansion process. 

2. Equations (Bl) and (B2) are solved simultaneously to give yp and tp, the 
coordinates of point P. 

3. The coordinates of point D must now be established to indicate where a particle 
path through point P intersects line CA or CB. This step is necessary because the rate 
equations and energy equations must be evaluated along the line DP. Equation (B5) gives 
the particle path as yptp = Yj^d* Three possible locations of point D must be consid- 
ered. If yptp > y c t c , point D is on line CB. Since for the frozen expansion CB is a 
straight line, it seems reasonable to treat CB as a straight line over the small interval 
considered herein. Writing the equation of the line CB and solving this equation simul- 
taneously with the condition yptp = Yj^b yield the values of the desired coordinates yj) 
and t D . A linear interpolation between points C and B is used to obtain values of the 
properties, compositions, and velocity at point D. If yptp < y^t^, point D is on line CA. 
The line CA is of the form, yt^ = Constant for the frozen expansion, so this form 
appears reasonable for use herein. Thus, by assuming Yj^b^ = t^Yc^C^ + Ya^A^) anc * 

by using the condition yptp = y^D* the coordinates yj) and tj) are readily found. A 
linear interpolation between properties at points A and C is used to evaluate the flow 
quantities at point D. If yptp = Yb^D> points D and C are identical and properties at D 
are defined. 

4. The rate expressions for the rate processes considered are evaluated at points A, 

B, P, and D for use in the compatibility and rate equations. Thus and can be 

calculated for all species by using equations (35) and (37). 

5. The compatibility equations, equations (B3) and (B4), are solved simultaneously 
to obtain improved values of pp and Up. The value of tp is determined by step 1. 

6. Rate equations (B7) and (B8) are applied along the particle path from point D to 
point P to obtain improved values of xi and cq at point P. 

7. The finite-difference form of the energy equation, equation (B6), is applied along 
particle path DP to yield the following new value of hp: 



(B9) 
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8. A new temperature at point P, found by using the definition of enthalpy, is 


n 



(BIO) 


9. The thermal equation of state is utilized to calculate an improved density at 
point P. 

10. Improved values are thus available for the thermodynamic properties, composi- 
tions, and velocity at point P. The improved values are used with the previously assumed 
values in a second-order correction scheme to obtain an improved approximation of the 
values at point P. The calculation returns to step 2 and the procedure is repeated with 
these new values. The iteration process continues until successive calculated values are 
within a prescribed degree of agreement. Inasmuch as temperature is more sensitive to 
nonequilibrium effects than other properties, successive temperature values are com- 
pared to establish the convergence of the iteration. When successive values are within a 
prescribed range of agreement, the calculation proceeds to the next grid point. 


Calculation of Near -Steady Region 

A calculation procedure has been established for computing the flow in the expansion 
fan. The near-steady region, the region between the tail of the expansion fan and the 
piston face, may also be calculated by using the same procedure with slight modifications. 
Two significant changes are necessary. First, it is convenient to use the b,t coordinate 
system instead of the y,t coordinate system. This change affects only the equations 
specifying the characteristic directions. Second, the boundary condition of constant- 
velocity flow is imposed at grid points on the piston face. At these points, up is fixed 
and equation (B4) is adequate for obtaining an improved value of Pp; thus step 5 of the 
calculation procedure is unnecessary. Interior points are calculated by using essenti- 
ally the same procedure as that used for the expansion fan. 
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TABLE L- CHEMICAL RATE CONSTANTS 


(a) Dissociation reactions 


Constant 

Value for reaction — 

0 2 + O - 30 

0 2 + 0 2 — 0 2 + 20 

0 2 + N 2 - N 2 + 20 

Dj, cm 2 / mole-sec 
Ej, cm 2 / mole-sec 
Fj, cm 2 / mole-sec 

0.905 x 1020 
-1.0 
59 395 

0.3258 x 10 20 
-1.0 
59 395 

0.724 X 10 19 
-1.0 
59 395 


(b) Recombination reactions 


Constant 

Value for reaction - 

30 — 0 2 + O 

0 2 + 20 - 20 2 

N2 + 20 — O2 + N2 

0.60334 X 10 16 
-0.5 
0 

Dj, cm6/mole 2 -sec 
Ej, cm6/mole 2 -sec 
Fj, cm6/mole2-sec 

0.75417 X 10 17 
-0.5 
0 

0.2715 X 10 17 
-0.5 
0 


TABLE H. - VIBRATIONAL CHARACTERISTIC TEMPERATURE 
AND RATE CONSTANTS 


Constant 

Value for species — 

°2 

n 2 

01, °K 

2270 

3390 

Ai, atm-sec 

0.1424 

560.61 

Bi 

0.3718 

-0.0995 

Ci 

179.33 

163.63 


32 












y = a 

y-t plane J o 


Figure 2.- Boundary conditions in the three coordinate systems. 
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Dimensionless velocity, 



Figure 8.- Velocity variation along specific particle paths in an unsteady expansion. 



Dimensionless time. 



Dimensionless modified Lagrangian coordinate, y 

Figure 10.- Complete grid network. (M.N) denotes a general point. 
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